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An Integral View of Fast Shocks around Supernova 1006 
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Supernova remnants are among the most spectacular ex- 
amples of astrophysical pistons in our cosmic neighbor- 
hood. The gas expelled by the supernova explosion is 
launched with velocities ~ 1000 km s ' into the ambient, 
tenuous interstellar medium, producing shocks that excite 
hydrogen lines. We have used an optical integral-field spec- 
trograph to obtain high-resolution spatial-spectral maps 
that allow us to study in detail the shocks in the northwest- 
ern rim of supernova 1006. The two-component Ha line 
is detected at 133 sky locations. Variations in the broad 
Une widths and the broad-to-narrow line intensity ratios 
across tens of atomic mean free paths suggest the presence 
of suprathermal protons, the potential seed particles for 
generating high-energy cosmic-rays. 

Supernova remnants, the expanding shells of material 
created in a stellar explosion, are astrophysical laboratories 
for studying non-thermal physics and high-velocity shocks, 
and have been scrutinized over a broad range of wavelengths. 
The signatures of non-thermal electrons are typically man- 
ifested in the X-ray, y-ray and radio range of wavelengths 
111]. Complementary to these observations, hydrogen emission 
from "Balmer-dominated shocks" ||2l|3l|4|], around supernova 
remnants of thermonuclear origin, directly probe the proton 
populations Q]. Until now, all studies of Balmer-dominated 
shocks have typically used conventional spectrographs that 
yield rich spectral information but limited spatial information. 
Here we report on Balmer-dominated shocks around a super- 
nova remnant using integral-field unit (IFU) spectroscopy, a 
technique that produces a three-dimensional "data cube": two 
dimensions of space (across the sky) and a spectral dimension. 
We selected supernova (SN) 1006 as our target, because it has 
a long history of serving as a laboratory for stud^ng non- 
thermal physics and high- velocity shocks laSHIilGSIlillSl- 

Long-slit spectroscopy, utilized in previous studies of 
Balmer-dominated shocks f\3i, cumulatively measures the Ha 
emission emanating from regions much larger than the charac- 
teristic length scale: the mean free path for interactions between 



hydrogen atoms and electrons or ions, L^fp ~ l/ncTce (with « 
being the pre-shock number density and cr^e ~ 10"'^ cm' de- 
noting the typical cross section for charge exchange). Using 
the inferred range of densities for the northwestern rim of S N 
1006 of n from 0.15 to 0.40 cm^^ HI El Ell, we estimated 
that the 0."67 size of a pixel in our IFU observations corre- 
sponds to ~5Lmfp. The high spatial resolution allowed us to 
separate out the contributions of bulk motion versus thermal 
velocity to the measured line widths, an issue which has lim- 
ited the interpretation of previous observations. It also enabled 
us to study the spatial variation of the line widths and ratios 
and hence changes in microphysics of a shock across several 
atomic mean free paths. To buttress the second point, we have 
intentionally chosen a field of view that zooms in on a region 
of simple geometry, implying that any substantial spatial vari- 
ation detected cannot be due to variations in density caused by 
geometric or projection effects. 

Figure [T] illustrates the general configuration of our ob- 
servations performed with VIMOS (Visible Multi-Object 
Spectrograph) in the IFU mode on the VLT. We have focused 
on the northwestern rim, because it produces the brightest emis- 
sion from a Balmer-dominated shock. As a first analysis of 
our data, we divided the shock structure into four strips and 
bin up the data along each strip. The two-component Ha line 
is convincingly detected. To first order, the broad line width 
(FWHM) W yields the shock velocity v.,, whereas the level of 
energy equilibration between electrons and protons in the post- 
shock gas introduces small corrections 11171 llSll . Because the 
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broad line profile is a direct probe of the velocity distribution of 
the post-shock protons, and the narrow hne profile traces pre- 
shock hydrogen atoms, the ratio of broad to narrow line intensi- 
ties lb I In contains information on how energy is shared between 
the pre- and post-shock regions. If no energy is shared, mod- 
els which include the basic shock and atomic physics correctly 
translate the measured W and hi In values into the inferred 
and p values 111 El Eol 12111 . where p is the ratio of elec- 
tron to proton temperatures in the post-shock gas. Unusually 
low hi In values (< 0.7) indicate that suprathermal particles 
from the post-shock region are traveling upstream into the pre- 
shock gas and depositing energy via atomic interactions (excita- 
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Figure 1 : VIMOS-IFU spectroscopy of the shock front in the remnant of SN 1006. The top-left panel shows a composite image of the full remnant (a30' in diameter), 
combining data from the Very Large Array and Green Bank Telescope (red; NRAO/AUI/NSF/GBT/VLA/Dyer, Maddalena & Comwell), Chandra X-ray Observatory 
(blue; NASA/CXC/Rutgers/G. Cassam-Chenai, J. Hughes et al.), 0.9 m Curtis Schmidt opdcal telescope (yellow; NOAO/AURA/NSF/CTIO/Middlebury College/F. 
Winkler), and Digitized Sky Survey (orange and light blue stars). The green box indicates the region covered by the Hubble Space Telescope (HST) Ha nan'ow-band 
image shown in the top-middle, with subsequently the magenta box indicating the region observed with the VIMOS-IFU. The top-right panel shows the reduced data 
cube collapsed in wavelength around the Ha-line, recovering the shock front. The crosses with four different colors indicate the pixels for which the spectra have 
been combined to produce the spectra shown in the four panels at the bottom. In each of the bottom panels, the dashed black lines show the best-fit double-Gaussian, 
with parameters given in the legenda: the intensity ratio of the broad to narrow component It 1 1,,, the velocity offset between the broad and narrow line centroids 
AV (in kms"'), and width of the broad component W (in kms"'). The reduced values above unity along with the differences between the observed spectra and 
their best-fits (shown below the spectra with an olfset of -0.3 for claiity), indicate that non-Gaussianity is present. Most of the reduced values above unity come 
from the mismatching near the line core, as can be seen from the zoomed-in region (-1000,1000) kms"' in the top-right corner of bottom panels. On the horizontal 
axis is shown only the fitted region of the spectra, while the y axis shows the flux in units of 10"'* erg s"' cm"^ A"' rescaled with the respect to the blue panel by 
factors 2 and 0.5 for the yellow and red panel, respectively. 



tion, ionization and charge exchange), thus acting as precursors 
1 22l l23l 12411 . The loss of energy decreases the shoclc velocity 



and hence W. The increased flux of particles into the pre-shock 
region leads to enhanced excitation of the hydrogen atoms via 
collisions and increases /„. The binned regions in Figure [Hal- 
ready hint at this phenomenon, because hi In ~ 0.6-0.7 in two 
of the strips (yellow and red). The broad line exhibits subtle 
deviations from a Maxwellian profile in the line core, further 
supporting the presence of suprathermal particles in the shock. 

We refined our approach by using the technique of Voronoi 
binning |25]. We combined neighbouring spectra into spatial 
Voronoi bins until S/N of about 40 was reached (see lower-right 
panel of Figure |2]i. We fitted double-Gaussians to the detected 
Ha lines as shown in Section 2 of the Supplementary Materials 
(hereafter SM) and derived W and /;,//„ values, thus producing 
the maps in the top-left and top-middle panel of Figure |2] In 



addition to the lack of complex density variations, we demon- 
strated that the viewing geometry of the shock is simple: by 
measuring the velocity shift AV between the centroids of the 
narrow and broad Ha line components, we derived the viewing 
angle via the relation AV - Sv, cos 0/4. Our measurements 
indicate that Q k 80°-90°. Therefore, the shocks we observe 
are mostly edge-on as expected, but the rim seems slightly S- 
shaped along the line of sight as suggested by il4ll . 

Surprisingly, we detected significant spatial variations in W 
and //,//„ across length scales 10"~70Lmfp, despite the simple 
geometry of the shocks, suggesting that they arise from varia- 
tions in the microphysics rather than density. Within the bright 
rim the variations in W are of order 10-20%, significantly larger 
than the individual measurement uncertainties (see Section 2 
and Table S 1 in SM). The density variations of 20-40% required 
to explain the detected variations in W are much larger than 
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Figure 2: Two-dimensional spatial- spectral maps of various properties associated with the shock front in the remnant of SN1006. The bottom-right panel shows 
the spatial Voronoi binning to reach a minimum signal-to-noise S/N ^ 40. The maps in the top row from left to the right: the broad line width W (in km s" ' ), the 
broad-to-nan'ow line intensity ratio It/I,,, and the shock inclination angle 6 (in degrees). Shock velocities i', and electron-to-proton temperature ratios /? shown in 
the bottom-left and bottom-middle panels, respectively. 



expected on these scales, and also incompatible with the un- 
changed smoothness of the shock over two decades of imaging 
observations ll9[ [l4ll . 

The low lb /In values are found at all distances from the inner 
rim, reaching as low as 0.4 at some locations (Figure |3}. We 
used the models without non-thermal physics to infer the val- 
ues of Vj and f3 lfl9ll . Most (about 85%) of the binned data are 
not accounted for by the model because of the low /;,//„ values, 
thus motivating the need for models which include suprather- 
mal particles and cosmic rays ll26ll . Combining the derived 
shock velocities with the proper motion measurement from 
we obtained estimates of the heliocentric distance to SN 1006. 
Clearly, there has to be a unique distance, likely given by the 
upper points around ~ 2 kpc (see also Section 3 in SM). Most 
points, however, are underestimating the distance due to loss 
of energy from the broad line component resulting in too low 
inferred shock velocities. 

These low values and variations of W and Ii,/I„, as well as the 
potential non-Gaussian contributions to the broad Ha lines, de- 
mand an explanation. We examined three possibilities. The first 
is that "broad neutrals" are acting as a precursor These 
are secondary populations of "hot" hydrogen atoms which are 
produced when post-shock protons capture an electron from a 
pre-shock hydrogen atom — the subsequent excitation of these 
broad neutrals produces the broad Ha line component detected 
in our observations ||2> El B E^l- Broad neutrals will warm 
the pre-shock gas and produce a third Ha line component of 



intermediate width, which may account for the narrow-plus- 
broad-line double-Gaussian being unable to fit the data within 
the measured uncertainties (i.e., not reaching a reduced x~ of 
unity). The strongest argument against broad neutrals is that 
they can only act as a precursor over an atomic mean free path, 
which is at odds with our observational result that the variations 
and low values of W and lb /In extend over distances » Lmfp- 

The second possibility is that pre-shock hydrogen atoms 
may cross the shock front, eventually become ionized and be- 
come protons gyrating along an ambient magnetic field line. It 
has been shown that these "pick-up protons" settle into a bi- 
spherical distribution which introduces a non-Gaussian contri- 
bution to the broad Ha Une core ll29[l . This explanation re- 
quires the magnetic field to be ordered on small length scales 
(<& Lmfp), whereas turbulent magnetic fields will result in a 
broad Ha line that is approximately Gaussian. 

Both explanations are disfavored at low pre-shock neutral 
fractions, which is the situation in SN 1006 with a pre-shock 
neutral fraction of about 0.1 lll3ll . This is because both broad 
neutrals and pick-up ions require the pre-shock gas to contain a 
substantial population of hydrogen atoms (relative to electrons 
and protons) in order to initiate the process. In the limit of a 
fully ionized pre-shock gas, no broad neutrals or pick-up ions 
may be produced. 

The explanation we favor is that the post-shock proton pop- 
ulation includes a non-thermal sub-population of protons — 
suprathermal protons (that are not pick-up protons). Such an 
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Figure 3: The three panels show all measured (with error bars ) broad line widths 
W, broad-to-narrow intensity ratios Ib/I„, and heliocentric distances D from 
combining the proper motion measurement with the shock velocities. Data are 
ordered in increasing distance from the inner lim (black/purple) to the outer 
rim (orange/red color), shown in arcsec at the top and in units of lO'* cm at the 
bottom horizontal axis. The inner rim coincides with the inner edge of the blue 
rim in the top-right panel in Figure [T] The dashed horizontal lines indicate the 
measured IV=2465.76 km s"' and /fc//„=0.79 values from collapsing all spectra 
of the pixels on the shock front, and from there the inferred D= 1 .72 kpc. Data 
points marked with a filled circle are those for which a valid model solution is 
found. The triangles (squares) indicate that due to low (high) It /In values the 
models hit the lower (upper) boundary limit in /J of 0.01 (1.0). 

explanation requires no assumption on the magnetic field ge- 
ometry or pre-shock neutral fraction. The lack of non-thermal 
X-ray and TeV y-ray emission in the northwestern rim of the 
remnant indicates that if there is a CR precursor, the injection 
of electrons is little and the particle energies are low (< IMeV), 
but CR acceleration can still happen. 

Our pilot project demonstrates the feasibility of using 
integral-field spectroscopy to observe and study the micro- 
physics of high-velocity shocks around supernova remnants. 
The resulting high spatial resolution mapping of the Balmer 
dominated shocks in the northwestern rim of SN 1006 suggests 
the presence of suprathermal protons of energies 10-100 keV 
which can seed high-energy cosmic rays. 
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Supplementary Materials 

This Supplemental Material consists of three parts. The first 
part describes our intergral-field spectroscopic observations and data 
reduction. This is followed by an assessment of the accuracy in fitting 
the Ho- emission line profile to extract the shock properties. Finally, 
we discuss estimates of the heliocentric distance to SN 1006. Figures 
S 1 through S3 support the Ho'-profile fitting, with resulting parameters 
and uncertainties as well as shock properties given in Table S 1 . 

1. Integral-field Spectroscopic Observations and Data Re- 
duction 

Our observations of the northwestern rim of SN 1006, were per- 
formed in queue mode on the nights of 2010 April 6, 11, 12, 13, 
and May 10, using the VIMOS (Visible Multi-Object Spectrograph) in 
the IFU (integral field unit) mode on the VLT. The IFU was placed 
approximately 15' away from the remnant's center, at coordinates 
0-72000 = 15''02'"13'.5, (5j2ooo = -4r45'22". The overall exposure 
time of 6 hours accumulated in 10 exposures, all taken at air mass 
< 1.60 and seeing < 1".60, enabled us to collect enough photons to 
reach the required signal-to-noise ratios (S/N). A dither pattern with 
small offsets between the individual exposures was applied in order to 
account for dead fibers. We have used a spatial scale of 0".67 per pixel 
in combination with the HR-Orange grism providing a field-of-view 
of 27"x27"as well as a resolution of Rsi2650 within the wavelength 
range of 5250 A < ^ < 7400 A. 

The data has been reduced with our own dedicated data reduction 
pipeline which was developed in particular to deal with the rather com- 
plex VIMOS instrument in IFU mode. It is based on several individ- 
ual scripts written in the software language python. VIMOS is made 
of four identical optical quadrants, each representing a completely in- 
dependent spectrograph. Each quadrant is at first reduced separately, 
including the following standard procedures: bias subtraction, stray- 
light correction, interactive fibre-identification and -tracing, spectra 
extraction and wavelength calibration using the arclamp frames. For 
the extraction of the fluxes the "optimal extraction algorithm" was ap- 
plied isoll . which was designed to deal with dense packed IFUs to cor- 
rect for the cross talk between adjacent fibres. The associated contin- 
uum lamp frames were used to construct corresponding fibre-flatfields 
which were used to normalise the wavelength-dependent throughput 
of each fibre. Additionally, VIMOS also suffers from substantial dis- 
tortions on the raw frames due to instrumental flexure. Thus, in or- 
der to guarantee the best possible extraction results the fiber traces 
were smoothed and corrected for shifts between science and calibra- 
tion frames. This included a 2nd order correction of the absolute posi- 
tions of the traces as well as an adjustment of the wavelength solution 
by using the available strong night-sky emission lines within the ob- 
served wavelength-range. The raw frames were cleaned for cosmic 
rays using the newly developed software package PyCosmics ifsill . 
which represents an optimized algorithm to detect their artifacts in IFU 
raw-data. 

After the basic reduction steps the four individual science frames 
for each observing block were combined and rearranged into three di- 
mensional datacubes, using the dedicated VIMOS lookup-tables. Flux 
calibration was achieved by reduction and extraction of photometric 
standard stars which were observed using the same setup as the science 
data. As the VIMOS-IFU does not have "sky-dedicated" fibers, it is nor- 
mal practice to do separate offset "sky exposures", but in this case the 
sky background could be reliably measured from the "empty" down- 
ward triangular part in the field-of-view (see top-right panel of Fig- 



5 

ure I). The last step involved the correction of the barycentric Doppler 
shift, i.e., the offset in velocity space caused by the movement of the 
Earth around the Sun. The observations were accomplished within a 
time-window of approximately one month, causing a velocity shift of 
about 20 km s"' which corresponds roughly to one pixel on the de- 
tector. In order to readjust this offset we measured the center of the 
narrow Ho- emission line and shifted it to the laboratory wavelength of 
6562.8 A. 

The Ha line is found in the middle of the HR-Orange grism wave- 
length range of 5250-7400 A. This spectral range is broad enough to 
both easily fit the broad component of the Ho- line with width of the 
order of 50 A, and to still leave plenty of spectral range for an accurate 
determination of the continuum. Finally, we used the precise header 
information of the telescope-pointings to compute the spatial-offsets 
between each observation and combined all ten individual datacubes 
within one final cube, containing in total 2162 spectra. The instrumen- 
tal resolution of about llOkms"' (FWHM) around the wavelength of 
Hcf is more than sufficient to measure the width of the broad line of 
2000-3000 kms"', as well as deviations from a Gaussian profile. The 
width of the narrow line, expected to be about 20 km s"', cannot be 
resolved, but the narrow line intensity can be accurately measured, es- 
pecially given that the narrow line has a significantly higher S/N than 
the broad line. 



2. Accuracy of the Ha-line fitting 

We fitted the Ha-line with two Gaussians, one for the narrow com- 
ponent and the other for the broad component, after convolving with a 
Gaussian instrumental profile. Hint on the non-Gaussianity around the 
line core of the broad component that we detected is a supporting, but 
not the main piece of evidence of suprathermal protons presence. 

In order to investigate possible non-physical contributions to the 
Ha-line profile, we compared the continuum emission profile across 
wavelength in the region downstream and upstream of the shock. As 
shown in Figure ISTI we combined spectra in the downstream and up- 
stream regions. The continuum has the same profile in these two re- 
gions, and does not affect the Ha-line profile. In these regions we did, 
however, detect diffuse narrow-line Ho- emission, with surface bright- 
ness, 10-30 times lower compared to the emission in the four regions 
in Figure I . 

Next, to verify that the instrumental profile is Gaussian we fitted 
sky lines that went through the same data reduction steps (except for 
sky subtraction). In Figure |S2] we show five sky lines around Ha line 
fitted with a Gaussian. The fits together with the differences between 
the observed and fitted sky lines show that the instrumental profile is 
indistinguishable from Gaussian. 

We have also checked if the residuals might be coming from the in- 
trinsic variations among the individual spectra that went into the com- 
bined spectra. The offsets between velocity centroids are much smaller 
than the "extra-core" size, and thus rule out this possibility. Also, the 
individual bins have veiy similar narrow line widths that are always 
unresolved, but much smaller than the "extra-core" size. 

Finally, the uncertainties in measuring the broad line width W and 
broad-to-narrow intensity ratio /;,//„ were determined in various inde- 
pendent ways, taking into account the varying noise per pixel in wave- 
length as given by the derived error spectrum per pixel. Two of the 
methods were very robust and led to consistent error estimates even 
in case of spectra with signal-to-noise significantly lower than our tar- 
get S/N X 40. In the first method, we created 1000 mock spectra 
per Voronoi bin through Monte Carlo sampling of the observed (com- 
bined) spectrum, and fitted each mock spectrum with a double Gaus- 
sian. The resulting distributions of values in W and It,/I„ shows that 
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both quantities are uncorrelated and their means are the same as from 
the fit to the observed spectrum. Their standard deviations, however, 
are slightly larger than those from the second method in which we var- 
ietQ W and 4 /l, on a grid around their best-fit values until changed 
by an amount corresponding to the 68% confidence level given the de- 
grees of freedom. This is expected, as the first method, apart from the 
added noise from the error spectrum, also includes noise already in- 
herent in the observed spectrum, increasing the reduced values of 
the mock spectra fits above unity. Re-scaling the reduced values to 
unity leads to uncertainties that are fully consistent with those inferred 
from the second method. The double-Gaussian fits to the Ha-line pro- 
file in all 133 sky positions are shown in Figure IS3l The resulting best 
fit parameters and uncertainties are given in Table S 1 , along with shock 
properties inferred using van Adelsberg models. 



3. Heliocentric distance to SN 1006 

When converting from angular to physical length scales, l"w 
2.58x10'* cm, we adopted a heliocentric distance of 1.72 kpc, derived 
as follows. We first summed up all pixels on the shock front, and ex- 
tracted the width of the broad component W and intensity ratio It/I,, 
by fitting the Ha line with double-Gaussians. The resulting broad line 
width and broad-to-narrow intensity ratio are W = 2466 kms"' and 
h/I,, = 0.79. Using the van Adelsberg et al. model, this yields a shock 
velocity of v.,=2288 kms"'. Combining this shock velocity with the 
proper motion measurements of 280 mas yr"' [9] yields a heliocentric 
distance of 1.72 kpc. 

Previously inferred distances to SN 1006 from the Ha observations 
in the northwestern rim, 2.18 kpc [9] and 1.6 kpc [19], were based 
on long-slit spectroscopic observations from [13]. The discrepancy in 
the inferred distances comes from different shock models applied. The 
latter one uses the same van Adelsberg et al. model, yielding a lower 
shock velocity for the same W = 2290 km s"' and 4//,, = 0.84, pri- 
marly due to the broad neutral velocities contribution to the relative 
speeds in fast neutral-ion interactions. 

Analysis of the optical spectrum of the Schweizer-Middleditch 
star 13211 which lies almost at the same line-of-sight as the center of 
SN 1006, sets the upper limit on the distance of SN 1006 at 2.1 kpc. 
The observed ejecta expansion at 7026 kms"' ifssll , and the require- 
ment that this material lie within the remnant places a lower limit of 
1.6 kpc. We actually expect the distance to be closer to ~ 2 kpc to 
explain the upper envelope of the points in Figure 3, but a model that 
properly takes into account the non-thermal physics is needed before 
we can robustly measure an intrinsic shock velocity and correspond- 
ingly obtain a secure estimate of the heliocentric distance. 
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Keeping the total intensity /„ + /;,, narrow line width and centroids of the 
nan'ow and broad component fixed to their best-fit values. 
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Figure S 1 : The left panel shows the reduced data cube collapsed in wavelength around the Ha-line with blue and red crosses indicating the pixels for which the 
combined non-sky-subtracted spectra are shown in the right panels. Comparing the continuum emission across wavelength, we confimi that the emission profile is 
the same in the upstream and downstream region of the shock. 
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Figure S2: The top panel shows a non-sky-subtracted spectrum with five sky lines around Ha, indicated by vertical dashed lines and numbered. Next, each row 
shows one sky line fitted with a Gaussian for the four different regions consistent with the ones in Figure 1 . The fits along with the differences between observed 
and fitted sky lines (with an offset downward for clarity) show that the instrumental profile is Gaussian. 
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Figure S3: Montage of the Ha line profiles in each of the 133 Voronoi bins. The horizontal axis represents wavelengths and vertical axis is in the relative units. 
Over-plotted are the double-Gaussian fits to the narrow and broad line components. No smoothing has been applied to the data. 
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Table SI: Properties of the shock front in the remnant of SN 1006 for 133 spatially binned locations. Columns 1-5: number of the (Voronoi) bin, x and y coordinates 
of the bin centroid, number of combined pixels, and signal-to-noise ratio. Columns 6-10: measured values of the broad component W and the broad-to-narrow line 
intensity ratio It/I,, (both with estimated uncertainties based on;^'^ 68% confidence levels), shock inclination angle (with typical uncertainty of 1.6° coming from 
a typical error of 47 kms^'in AV), and the reduced of the double-Gaussian fit to the observed Hty-line. Columns 10-11: shock velocity v, and electron-to-proton 
ratio fS (in base-10 logarithm) of the best-fit shock model; a long dash indicates when the shock model hits the adopted boundary limits of log/? = (-2,0), Hkely 
because CR physics is missing. 
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